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SUMMARY 


In this report we describe the two-dimensional Lagrangian, incompress- 
ible Cartesian code, SPLISH, and the changes made to convert it for the study 
of flows in and around fuel droplets. The Lagrangian technique used in this 
study incorporates a general restructuring triangular mesh, which allows 
reconnection of vertices to eliminate grid distortions without adding 
numerical diffusion, "fliis technique is accurate at material interfaces even 
though the interfaces undergo convolutions and may evolve into multi- 
connected surfaces. 

New algorithms for surface tension and viscosity have been added to the 
basic fluid dynamics code. Surface tension is included as a jump in pressure 
across an interface by casting the surface tension forces in the form of a 
gradient of a potential. The surface tension algorithm is benchmarked by 
studying the oscillatory behavior of an n = ? normal mode. The viscosity 
algorithm for a general mesh is presented and tested by calculating the 
spreading of a viscous shear layer. 

We use the code to calculate the internal and external flows of 
oscillating and deforming kerosene droplets in an air jet. The surrounding 
air jet is initialized to laminar flow about a round kerosene droplet. 'Kie 
evolution of the droplet and jet are calculated from first principles, 
eliminating approximations for effective droplet size, wake effects or 
recirculation patterns. Results of the air-kerosene calculations are 
illustrated by sequences of frames from a computer generated movie of fluid 
particle positions. Both internal and external flows are shown as well as 
droplet distortion due to the relative flow. The algorithms needed to extend 
these calculations to compressible flows in three dimensions are discussed. 



INTRODUCTION 


Droplet combustion is a complex transient problem in multiphase flow. 
Particularly severe physical and mathematical approximations must be made to 
describe the detailed interactions between droplets and the external flow 
field in spray combustion models (Williams, 1973; Faeth, 1977,1983). For 
example, equivalent spheres are used to approximate droplet deformations, and 
empirical expressions are used to account for drag and convection. The 
effects of droplet breakup are included by using estimated breakup times and 
drop sizes after breakup. Quasi-steady flow approximations are used, and 
changes in the flow field due to droplet deformations, wake effects and 
droplet distortions due to the flow field are neglected. Finally, in most 
models the droplet concentration is assumed to be dilute since little is 
known about droplet-droplet or droplet-wake interactions. 

The need for these approximations arises directly from the difficulty in 
following several physically distinct regions as they interact with the 
external flow field, distort and separate or merge, A Lagrangian technique 
is well suited to accurately modelling the transport of these various regions 
since it easily and naturally calculates the advection of boundaries. 

Because the various regions, may severely distort and separate, the 
calculational grid must be able to self-consistently adapt to the physical 
flow. For these reasons the numerical technique used in this study is 
transient hydrodynamic modelling using a Lagrangian mesh. 

The calculations are performed using the fully conservative, two- 
dimensional Lagrangian finite difference method developed by Fritts and Boris 
(1979) specifically to handle multi -phase flow. The method is based on a 
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dynamically restructuring Cartesian triangular grid. Triangle sides are 
aligned on material interfaces. Since vertex movement is Lagrangian, the 
interface sides accurately track the movement of the interface due to 
advection. A triangular grid avoids the problems of mesh tangling 
encountered in Lagrangian methods using a quadrilateral mesh: individual mesh 
points are continually reconnected to account for the migration of fluid 
elements in the flow field. Since the number of grid lines meeting at a 
vertex is variable, the resolution can be altered non-dif fusively where 
needed (e.g., around a region of droplet distortion) without affecting the 
resolution in other areas of the computation. This is a major step forward 
in the computation of droplet flows because the Lagrangian technique allows 
for the evolution of the grid to multiply-connected regions. Thus, according 
to the flow conditions, droplets can break up and shatter. 

Previously, the Lagrangian restructuring triangular grid technique has 
been applied to a number of incompressible fluid flow problems including 
calculations of nonlinear waves, flows over solid obstacles, Kelvin-Helmholtz 
and Rayleigh-Taylor instabilities, Couette flows and Taylor vortex flows 
(Fritts, 1976,1976a, Fritts et. al., 1980,1981, Emery et. al., 1981). This 
report details the research performed in adapting this technique to the 
study of droplet flows and presents calculations of the flow about droplets 
from first principles, without resort to models for droplet distortion, 
breakup, wake effects, drag or convection. The goals of the project are to 
extend this work further to build a comprehensive model for droplet 
combustion. This paper is a final report detailing the form of the model and 
indicating the results obtained in testing the model on purely hydrodynamic 
flows about droplets. Although the effects of combustion are not included in 
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this report, the model was constructed in such a way as to be able to account 
for the additional complexity introduced by combustion processes. 

Consider the burning of an isolated heated droplet. The droplet heats up 
from its surface inward. Depending on the temperature, fuel, and other 
ambient conditions, the droplet may develop substantial internal gradients. 

If there is enough convection, recirculating flow develops within the 
droplet. As the surface evaporates, the locus of fuel and oxidizer at the 
fuel rich limit expands radially outward from the droplet. After the correct 
chemical induction time, the mixture may ignite outside this locus. For 
burning to continue, heat from the oxidizing gas and products must diffuse in 
past the products to continue the evaporation process, and fuel must continue 
to diffuse outward toward the flame region. Outside of the flame region the 
expanding oxidizing gas may sweep away unburned fuel vapor as well as 
combustion products. The Lagrangian technique improves upon phenomenological 
models because the advective motion within these several distinct regions are 
accurately calculated without prior assumptions of flow patterns and without 
approximating transient flows by steady-state or quasi steady -state flows. 
Adjustments of interface positions due to thermal, mass or energy fluxes 
across a triangle side may be calculated in a separate step using a 
conservative integral technique. 

Although the physical processes described above are complex and highly 
nonlinear, the description is still idealized. A practical combustor must 
establish a relative velocity, V^j, between the fuel droplets and the hot, 
oxidizing gas. Unless becomes substantial, surface tension keeps the 
droplet essentially spherical as it shrinks due to evaporation. For large 
Vjj, the droplet distorts and may even shatter. As the droplets move, a 
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boundary layer of vaporized fuel and hot -gas develops which creates drag, 
decreases and may introduce circulation within the droplet. Recircu- 
lation patterns may also develop outside the droplet, influencing both the 
boundary layer and external flows. The use of models to account for the 
droplet deformations and viscous effects may be avoided by incorporating 
algorithms for surface tension and viscosity into the hydrodynamics. Ihis - 
report describes these algorithms and the benchmarks used to validate their 
accuracy. 

The shape of the flame surrounding the droplet depends strongly on the 

/ 

distorted external velocity field. As the Reynolds number increases, the 
flame shape around the droplet changes from an envelope to a wake flame. The 
site at which energy is released therefore changes, and this in turn, 
readjusts the external flow field, local evaporation rates and species 
concentrations. To this picture must be added the interaction of burning and 
evaporating droplets. For close droplet spacings, the fuel-rich limit may ~ 
extend around all the individual droplets and the resulting flame strongly 
resembles a diffusion flame. Large droplets may still migrate through the 
boundaries of the flame, and the combustion characteristics of these droplets 
may strongly affect the concentration of combustion products. Finally, the 
physical boundaries of the combustor may be important, altering flow 
patterns, temperature gradients, or other important properties of the 
system. The Lagrangian technique is well suited to tracking the individual 
or merging droplets, the transient flame shape and the flow by irregular 
combustor boundaries. Future additions to the technique will allow for 
compressible effects, heat release, evaporation and chemical reactions (Oran 
and Boris, 1981). 
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An adaptive grid technique can be implemented in two distinct ways: 
recomputing the entire grid at each timestep or locally restructuring the 
grid to eliminate distortions. Either method requires the storage of 
bookkeeping information necessary to compute the finite-difference 
templates. However, the mesh distortions which develop during a single 
timestep are confined to a fairly small number of triangles, typically at 
most 5 percent of the grid, so that the grid restructuring method can be an 
order of magnitude more efficient. This remains true even on vector 
machines, despite the fact that grid restructuring is inherently a scalar 
computation. The computer code used in this report was implemented on a TI- 
ASC parallel processor computer using vectorized scans to test for locales 
where grid restructuring might be necessary. Scalar routines then performed 
the actual calculations for grid restructuring in those regions only. The 
hydrodynamics routines were all vectorized where possible. Efficient coding 
for the ASC was obtained through Fortran-callable assembly language routines 
which optimized the scalar fetch and store operations. The main drawback of 
the technique is that the resultant programs are machine specific and must be 
recoded to run on other machines. Because the programs are not transport- 
able, the code is not available for general distribution. Persons interested 
in obtaining program information or using the code should contact the authors 
directly. 

The complete computer code is roughly 13,000 lines long, including 
documentation. Execution times on the ASC are typically about .01 seconds 
per timestep per grid point, including program diagnostics and output in the 
form of two three-color movie films, individual frames on fiche, and fiche 
listings. A movie supplement of various film sequences discussed in the text 
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is available for use with this report. Movies are generated through an 
optimized package which outputs to a Dicomed film recorder, Tektronix 
terminals or Calcomp plotters. 

The calculations illustrated in the movie sequences and in the individual 
frames in the report all illustrate grids for a single droplet. Because the 
boundary conditions for all the calculations are periodic at the sides of the 
computational region and reflective at the top and bottom, the calculations 
represent an infinite series of droplets. However, most of the calculations 
made for this report terminate when the wake of the preceeding droplet 
impinges on the droplet following. This permits initial calculations of 
nearly isolated droplets. For the 125 micron drop size and the flot; speeds 
used in the calculations, droplet distortions and possible shattering were 
expected, and single droplet simulations of this behavior could be compared 
more easily to experimental observations. The gridding routines can generate 
initial grids for any desired drop size within an arbitrarily large mesh. 
Large differences in resolution are permitted so that the droplet interface 
and interior can be well resolved despite larger grid sizes in some regions 
of the external flow field. The current code is therefore capable of 
simulations of arbitrarily large or small droplets. 
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LAGRANGIAN HYDRODYNAMICS ON A TRIANGULAR GRID 


1 . General Approach 

In principle, a Lagrangian formulation of the hydrodynamic equations is 
particularly attractive for droplet combustion calculations. Each fluid 
element is tracked as it evolves through the. interaction with its changing 
environment and with external forces. Heat release, contaminant reactions 
and soot formation can all be represented locally, without resort to global 
models and without nonphysical diffusion. Conservation laws are simple to 
express since there are no fluxes out of the fluid element boundaries and the 
paths of the fluid elements themselves provide flow visualization. However, 
in all but the simplest flows the individual fluid elements deform, and these 
deformations are a severe hindrance to actually using a Lagrangian method. 

In numerical calculations, fluid element distortion appears as 
stretching, shearing and eventual tangling of the computational grid. 

Although the use of a general-connectivity triangular mesh eliminates 
tangling, the accuracy of a calculation may still deteriorate when there are 
abrupt local changes in resolution and when the high-order effects of 
deformations are not represented. Therefore it is very important to pay 
close attention to how well conservation laws are satisfied. For example, 
the accuracy of the finite-difference algorithms for a general mesh may not 
be sufficient to conserve quantities advected with the fluid elements if some 
flux is allowed to flow out of elements to maintain straight lines in the 
computational grid. In the following section we show hew exact conservation 
may be maintained. 
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The divergence and curl of the velocity field prescribe the kinetics of 
the field by specifying the local rate of expansion of the fluid, d, and 
local vorticity, C/ by 

V • V = d 


V X V = 1 • 

For incompressible flow, d = 0, and for irrotational flow, C = 0. 

For incompressible and irrotational flow in two dimensions, the velocity 
field is specified by a velocity potential 4’ and stream function i}»: 

_ li _ 

X 3x 9y 

= 

y 3y ~ 3x . 


( 1 ) 


( 2 ) 


These equations automatically satisfy the conservation of vorticity and mass, 
since 

• '1/ 

V •v=V«(Vx4j) = 0 

and (3) 

Vxv = VxV(fiE0. 

We would define finite-difference operators for divergence, curl and gradient 

which have these identical properties, and this requirement restricts the 

placement of variables. In particular, if 4* and 4> are to be assigned to the 

Lagrangian verticies, the velocities v must be specified at the 

centroids of triangles or the midpoints of line segments. Therefore the 

Lagrangian vertex velocities must be obtained by local averages. 

For example, the first of Eqs. 3 will be recast in finite difference 

notation. 1310 notation J. , , is the sum over vertices i around a 

^i(c) 

central vertex c. In such sums the sequence of vertices is assumed to be 

counter-clockwise around the central vertex. The quantity A. 

1 + 1/2 
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represents the area of the triangle having vertices (c,i,i+ 1 ). In Figure 1 
the area of triangle j is given by 


2», = (r3 - r^) X Ir, - r^) • z, 


( 4 ) 


where z is the direction out of the page. Similarly, for a scalar function f 
specified at each vertex and assumed piecewise linear within each triangle, 
the vector gradient of f (constant throughout the triangle j and 
discontinuous at the triangle sides) is given by 


zx(r -r ) zx{r -r ) zx{? -? ) 

2I * " *3— ST-^ 


( 5 ) 


3 2 3 

With this placement of variables the dynamics of the flow, as well as the 

kinematics, behave properly. That is, 

V*v = V»V(j)= V^({> (6) 

is a general triangular grid Poisson equation which may be used to solve for 

the local pressure. At the same time, the pressures generated by forcing all 

local divergences to zero cannot by themselves alter the local vorticities, 

due to Eg. ( 3 ). In finite-difference form Equation 6 becomes 

zx(f -r. ,) zx(f.-r ) zx(r. -r. ) (r. -r. ) 

n - V r, c 1+1 . X c , x +1 1 1 1+1 1 ^ 

1*1 2 A. 2 ». 2 ft. 1 " 2 

i(c) 1+1/2 1+1/2 1+1/2 

where A is the area of the vertex cell, defined as one third of the 
c 

sum of the areas of all triangles including that vertex. 

The accuracy of the numerical algorithms is determined by both the local 
resolution and connectivity of the grid. For the approach used here, the 
local connectivity of the grid and the resolution are both determined in part 
by the requirements that the matrix generated from the Poisson equation, 
Eq.( 7 ), remains diagonally dominant. With this restriction, convergence of 
an iterative solver for Eq. ( 7 ) is assured. The consequences of maintaining 
diagonal dominance will be given below in the discussion of grid 
restructuring algorithms. 
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2. Finite-Difference Algorithms 


In the previous section it was shown that as a result of specifying ip and 
4> on vertices, the velocities must be specified at triangle centroids. With 
this definition the vorticity C about any vertex is readily calculated. Two 
formulations of the basic incompressible hydrodynamics equations are 
accessible with these definitions. For a formulation the vorticity ^ is 
advanced at each cell for each timestep and the new stream function is 
obtained from a solution of Vx(Vx\|))= 5 . By Eg. (3) the divergence of the 
velocity field is identically zero. Alternatively, in a P-v formula- 
tion the changes in vorticity are zero by construction since VxVp=o. Then 
the new pressures are chosen to force the divergence of the velocity field to 
zero. Ihe P-v formulation has been the focus of previous work using 
this technique for several reasons. First, the formulation becomes more 
complicated in three dimensions. Second, boundary conditions around bubbles „ 
and cavities are more complicated in the version, particularly for 

droplet shattering and coalescence. Third, a self-consistent pressure field 
is usually desired even in the formulation, requiring an extra solution 
step. 


The P-v algorithm specifies pressures, velocities and positions 
at full timesteps. A split-step algorithm is used to integrate the 
velocities forward half a time step, advance the grid a full time step, and 
then advance the velocities the remaining half time step. 


-1/2 -o 6t , 6t ^ 

Vj - V. - --j-gy. 


( 8 ) 



-o 

V. 

1 


-n 

+ V. 
1 




(9) 


1 1 



-n 
r . 
1 


( 10 ) 


-o 
r . 
1 


+ 



t 





( 11 ) 


-n 

V 


D 


~1/2 6t ,n 6t 


gy. 


( 12 ) 


In these equations the subscript i denotes a vertex quantity, while the 
subscript j is used for triangle-centered quantities. Therefore the new 
vertex velocities appearing in the right-hand side of Eq. (9) are obtained 
from the area-weighted new triangle velocities found from Bq.(12) during the 
previous iteration. That is, Eqs. (8) and (12) advance the velocities 
according to the Lagrangian equations of motion: but since the grid is 
advanced at the half time step, a vertex velocity at the half time step must 
be found from the old and new triangle velocities. This implies an iteration 
over Eq.(9) through Eq. (12) to assure that the new velocities are indeed 
consistent with those used for the grid advancement in Bq.(IO). 

Equation (11) is the numerical expression of the change in the triangle 
velocities that must occur during the grid advancement if the vorticity is 
to remain contant for inviscid, homogeneous flow. This transformation is 
apparently a unique, but necessary, addition to Lagrangian methods to assure 
that vorticity is conserved. The actual form of the transformation has 
changed during this project, so further discussion will be deferred till 
later. 


12 



3, Adjusting and Restructuring the Mesh. 

The primary advantage of a restructuring mesh is the flexibility v/hich it 
permits for Lagrangian techniques in following long time solutions to 
complicated flows. Several types of local mesh adjustment and restructuring 
are used to maintain uniformity and accuracy of the discrete mesh 
representation. A mesh adjustment is a nonphysical movement or adjustment of 
the position of one or more vertices without changing the connectivity of 
mesh vertices. These adjustments are designed to regularize the mesh, and 
result in the effective transfer of fluid across triangle sides. 

Mesh restructuring, on the other hand, does not generally involve 
movement of vertices but generally a redefinition of the mesh connectivity. 
Simplification of the mesh under restructuring may also involve vertex 
addition and deletion, but the positions of all other vertices remain 
unchanged. Therefore adjustment and restructuring are somewhat orthogonal 
procedures, one leaving vertex positions unchanged and the other leaving the 
mesh connectivity unchanged. Since restructuring always involves the changed 
position of a triangle side, it can also incorporate the nonphysical flow of 
fluid across triangle sides. 

Both adjustment and restructuring represent departures from a purely 
Lagrangian description and threaten to introduce unwanted numerical diffusion 
into the method. To minimize diffusive and other errors, vertices and 
triangle sides lying on boundaries, surfaces and interfaces must be left 
undisturbed, and mass and moment^^m must be strictly conserved everywhere 
during both restructuring and adjustment. Although there are many schemes 
possible for mesh adjustment and restucturing, we have concentrated on a few 
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"primary" procedures from which more complex procedures may be developed. 
Since all conservation laws are satisfied for these simple procedures, all 
schemes built up from these primary procedures will also satisfy the same 
conservation laws. 

A triangular mesh can quickly become distorted through the migration of 
vertices in the fluid flow, particularly for shear flows. This situation is 
typified by regions of long, narrow triangles bordering more regular ones. 
Without restructuring this distorted mesh forces the computation of 
derivatives using vertices which are no longer the nearest neighbors, and 
quickly leads to inaccuracies, numerical instabilities and nonphysical 
behavior. Time-step errors also become severe because of the disparity in 
size of triangle sides. For extremely distorted triangles, triangle 
inversion becomes likely. Because of the severity of these problems, grid 

restructuring must be imposed continuously to insure the accuracy of the 

♦ 

numerical solution. 

On a triangular grid, every nonboundary line uniquely specifies its two 
bordering triangles. These triangles form a quadrilateral for which the 
included line is one of the two possible diagonals. Figure 2a illustrates a 
configuration for which the present diagonal (solid line) should be 
reconnected to the opposing diagonal (dashed line). One possible algorithm 
always chooses the shorter diagonal unless reconnection produces too large a 
disparity in triangle areas. This safeguards against reconnecting the 
diagonal of inverted quadrilaterals to produce a negative area triangle, as 
shown in Fig. 2b. 

The reconnection algorithm could instead be formulated to ensure diagonal 
dominance of the triangular grid Poisson equation (Bq.7) as mentioned 
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previously. Note that the coefficient of the <))^ term in Eq. 7 is 

I r . -r . I ^ 

^c ~ 4A. ' 

i(c) 1+1/2 

and is always negative. Ihe coefficient a^ of the (j>^ term is 


a . = - 
1 


(r -r. , ) *(r. -r. ) (r. -r )»(r.-r. ,) 

c 1+1 1+1 1 1-1 c 1 1-1 


"^i+1/2 


1 - 1/2 


Equation 14 reduces to 

where 6. _ and 6. ^ are the angles in the (i+1/2)th and 

i+1/2 1-1/2 ' 

(i-1/2)th triangles opposite the line from c to i. If the sum of 

and 6^ less than rr radians for each i, the matrix is diagonally 

dominant and normal iterative procedures for inverting Eg. 7 work well. If 


for any i 


6, , /_ + 6 > ir radians 

1 + 1/2 1 - 1/2 


then the line from c to i is reconnected to join (i+1 ) to (i,-1 ) . This 
procedure therefore chooses the other diagonal of the quadrilateral formed by 
the (i+1/2)th and (i-1/2)th triangles. Since the sum of angles in a 
quadrilateral is just 2 tt radians, then the new diagonal has opposite angles 
that sum to less than n radians. The new matrix coefficients generated by 
this connectivity again ensure that the matrix is diagonally dominant. Note 
that negative area triangles cannot form with an algorithm that requires that 
the sum of the opposing angles is greater than zero and less than ir radians. 
For all simulations presented in this paper, the algorithm enforcing diagonal 


dominance was used. 


The reconnection algorithm is complicated by the need to uphold 
conservation laws. To conserve vorticity locally, the new triangles defined 
by a reconnection have velocities constrained to those which leave the 
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vorticity about each vertex unchanged. The additional requirement that the 
momentum is locally conserved uniquely specifies the post-reconnection 
velocities for the two new triangles. The algorithm resulting from these 
constraints is reversible. Replacing the reconnected diagonal with the 
original diagonal redefines the initial two triangles with their identical 
original velocities. 

A further complication to the reconnection algorithm arises at material 
interfaces. Since triangle sides aligned along interfaces cannot be 
reconnected, diagonal dominance cannot be preserved for matrix coefficients 
from interface vertices by reconnection. Alternatively, a vertex may be 
added at the midpoint of the interface line so that the opposing angles are 
bisected by the lines drawn from the new vertex to the opposite vertices. 

This scheme assures diagonal dominance while increasing the resolution in the 
neighborhood of the interface. This algorithm is exercised when vertices 
close to the interface move toward the interface or when the interface 
becomes severely deformed. In either case, more resolution at the interface 
is generally required. 

The interface problem indicates that reconnection cannot solve all the 
mesh readjustment problems encountered in complex fluid flows. Two 
additional "primary" procedures are required: vertex addition and deletion. 

As discussed above, the addition of a vertex on an existing interface line is 
accompanied by the insertion of two new lines to form two new triangles. For 
a line on the boundary only one new line and triangle are added. The new 
vertex may be added anywhere along any interior, interface or boundary line, 
since later reconnections can be used to restore diagonal dominance. Two new 
triangle velocities must be specified, and these are selected in accordance 
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with the same conservation laws used for the grid restructuring algorithms. 

Vertices may also be added in the interior of any single triangle. 

Simple algorithms for vertex addition within a triangle tj^ically leave the 
grid motion unchanged. The three new triangles are circumscribed by the 
original triangle whose motion is usually constrained just as if the new 
vertex was not there, since the divergences, and subsequent pressure changes, 
are the same. To be effective, addition of a vertex within a triangle must 
be accompanied by the reconnection of at least one of the triangle sides. 
Physical variables centered at the vertex are found by interpolation just as 
in the case of addition on a line. 

Deletion of a vertex is performed by the inverses of these two processes. 
To delete a point within the interior of a subregion, reconnections are first 
made to isolate the point within a triangle. The vertex, three lines and two 
triangles are deleted, and the new physical variables are determined by 
averages over the old configuration. Subsequent reconnections enforce 
diagonal dominance. To delete a point on an interface, the interface must".-" 
first be realigned to its new, lower resolution position either by unphysical 
motion of the interface vertex or by changing the physical properties 
centered at one of the triangles. The inverse process can then be used to 
eliminate the vertex and redefine physical variables in accordance with the 
conservation laws. It should be noted that the reconnection algorithm may 
require the addition of a vertex at an interface for diagonal dominance, and 
may violate a resolution requirement stipulated in another part of the set of 
grid restructuring algorithms. Such conflicts can cause cycling through the 
addition and deletion algorithms unless the requirements are properly 
tailored to both requirements. 


17 



General grid restructuring algorithms are built from these primary 
functions, each of which is incorporated into a Fortran subroutine. The 
more general routines for grid restructuring are used to provide flexibility 
in setting resolution requirements and eliminating conflicts among the lower 
level functions. Input to these routines is in the form of a user specified 
resolution requirement limiting the maximum and minimum size for triangles 
and line segments and a CFL parameter which determines the time step from the 
flow speed. For the calculations presented in this paper these 
specifications were global in nature, although local specifications could be - 
used which would be based, for example, on material type, distance from 
interfaces or gradients in physical properties. Because the global 
resolution was determined by a range of acceptable sizes with the maximum 
near the initial values, finer resolution will persist wherever grid 
restructuring has occured, as is evident in the figures below illustrating 
computational grids for several different problems. Most of the finer 
resolution arises near interfaces, where the reconnection routines force the 
addition of vertices to ensure diagonal dominance and sufficient accuracy. 
Local resolution specifications were used only in the grid initialization 
algorithm, which is built from the same primary routines. 
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NUMERICAL ALGORITHMS AND RESULTS 


The basic two-dimensional hydrodynamics computer code was constructed in 
such a way that the finite-difference operators for divergence, curl and 
gradient exactly reflected the properties of the continuian operators. This 
construction assures conservation of vorticity and mass and provides a basis 
for determining the local grid connectivxty. The extensions to this code 
described below are all made in exactly the same spirit: the finite- 
difference approximations to both physical and dynamical processes conform to 
the continuum limit and conserve properly. 

The expansion of the basic triangular mesh code to droplet flows was 
programmed to occur in several stages. Each stage involved the development 
of new algorithms for particular additions to the physics being modelled or 
for necessary new numerical techniques, and the benchmarking of these new 
algorithms against relevant physical calculations. The following sections 
detail the progress achieved at the separate stages, the algorithms developed 
and the numerical results of the benchmark calculations. 

1 . Incompressible, Inviscid Flow about a Droplet without Surface Tension 

The first test problem for the code was a simulation of incompressible, 
inviscid flow about a cylindrical droplet with a density twice that of the 
background fluid. Gridding routines were written to position an arbitrarily 
large drop at the center of the computational grid for variable resolution 
inside and outside the droplet. Additional routines initialize the flow by a 
pressure pulse at the left boundary for the first half timestep or ramp up 
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the flow speed if the initial perturbation would cause adverse transients. 
Divergence calculations for advancing the pressure during the ramp-up time 
correct for the effect of the initial pressure pulse. For all later times 
the left and right boundary condition is periodic. Rigid wall boundary 
conditions are imposed at the top and bottom of the computational region for 
all times. The modeled problem is therefore a row of droplets in an 
impulsively started flow field. The droplet is gridded into 28 triangular 
computational cells in a total system of 552 cells, as seen in the first 
frame of Fig. 3. Only the triangle vertices are shown within the droplet. 

Figure 3 shows the triangular mesh at several times in the calculation. 
Pathlines of each of the vertices are plotted as a diagnostic, but are not 
included in this Figure. Early in the calculation a recirculation zone forms 
behind the droplet, compressing the droplet in the direction parallel to the 
flow. Flow within the droplet is initiated by this compression in a 
direction normal to the external flow, ihe bulges formed at the top and 
bottom of the distorted droplet are pulled around the recirculation zone by- 
the shear flow which is at a maximum at these points. The internal droplet 
flow is therefore driven by the compression set up between the front and rear 
stagnation points and by the high shear flow which extends around the top and 
bottom of the droplet and recirculaton zone. The interaction of the droplet 
back onto the external flow occurs primarily through the enlarged cross- 
sectional area presented by the droplet to the flow, which increases the size 
of the recirculation zone. Eventually, as seen in Figure 3, the droplet is 
squeezed into a thin layer coating the recirculation zone. The thinned film 
then shatters into several smaller pieces, first at the rear of the droplet 
and later in the more laminar flow toward the front of the droplet. 
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A Study of the pathlines of the Lagrangian particles shows that the flow 
is regular at all times despite the distorted shape of the droplet. 
Subtracting off the mean flow from the calculated flow field would show a 
large stationary recirculating double vortex. In a spherical geometry this 
recirculation zone would form a vortex ring, and the thinned droplet coating 
the ring would fragment in both the radial and azimuthal directions. The 
varying density of triangle vertices arises because higher resolution is 
required in the vicinity of the droplet interface. 

A second calculation was performed for a 10:1 ratio of droplet to 
external fluid density. The initial flow is quite similar and shows back 
flow at the rear stagnation point as well as internal droplet flow normal to 
the external flow. However, later droplet development is substantially 
altered, 'ftie more massive droplet is less easily deformed about the 
recirculation zone, and as a result the droplet grows more than the 2:1 case 
in the direction normal to the external flow. Tlierefore a more symmetrical '' 
f ront-to-back flow pattern develops. With no surface tension, there is no 
restoring force and no steady-state shape. The droplet grows normal to the 
flow until it is thinned sufficiently to break. Edge effects due to the 
proximity of the top and bottom of the computational region are clearly 
visible by the end of the calculation. 

The results of both tests agree qualitatively with existing theory and 
experiment. Because of the lack of surface tension, no quantitative 
comparisons could be made. The gridding algorithms were found to be 
sufficient to represent the droplet down to the desired resoluton as input to 
the calculation. For both calculations the grid adjusted itself 
automatically, i.e. without need for user intervention, to the changing 
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connectivity of the shattering droplet and to the flow about the stagnation 


points. Ihe Lagrangian pathline diagnostic was found to be effective in 
illustrating recirculating flow in a movie format. During the debugging 
phase of these calculations it was also necessary to develop a new diagnostic 
routine to zoom in and plot the grid and local field variables at points 
where the code indicated problems at gridding anomalies or in convergence of 
the Poisson solver. 

2. Surface Tension. 

The incorporation of surface tension into the code was a learning process 
in finite differencing over distorted grids. Surface tension is convention- 
ally cast into a finite-difference form by fitting vertices on the material 
interface to some parametric function from which an estimate of local 
curvature can be made. Once the curvature is known, a surface tension force 
is evaluated and used to accelerate interface vertices. This scheme fails 
for two reasons. First, the interface vertices are accelerated directly by 
surface tension forces evaluated on the vertices. Since velocities are 
centered on triangles in SPLISH, unless a secondary calculation is made, the 
velocity field sees the effect of the acceleration a half-timestep later, and 
as a result the pressure calculated within the droplet is inconsistent with 
that found from the surface tension formula. Secondly, since the pressure 
gradient forces and surface tension forces are not calculated in the same 
manner, niimerical error results which grows with each timestep. 

Both of the probleiiB mentioned above were eliminated by a new and unique 
formulation of surface tension in which a surface tension potential is used 
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to generate the forces. Since the pressure gradient forces are calculated in 
the same manner and on the same grid as those derived from the surface 
tension potential, exact balance can be achieved between the forces and 
static pressure drops across the interface agree exactly with theory. The 
surface tension force is then formulated as a gradient of a potential present 
only at the surfaces. Therefore both the "surface tension potential" and the 
pressure are dynamically similar, and the physical pressure drop across the 
interface must exactly cancel the surface tension forces. Since the surface 
tension is normal to the interface and opposes the pressure drop (Fritts, et. 
al., 1982), then the Vp x Vp terms which alter the vorticity are zero for the 
finite-difference algorithms. 

The finite-difference algorithms for surface tension are therefore quite 
simple in form. The surface tension forces are included through Laplace's 
formula for the pressure jump across an interface (Landau and Lifshitz, 1975) 


P - P = 0/R 
1 o 


(17) 


where p^ is the pressure just inside the droplet at the interface, P^ is 
the pressure just outside the droplet at the interface, a is the surface 
tension coefficient associated with the two media which define the interface, 
and R is the radius of curvature of the cylindrical droplet. The radius of 
curvature is positive at points on the interface where the droplet surface is 
convex (a spherical droplet is convex everywhere) and negative when the 
droplet surface is concave. Hiese pressure jumps are included in the Poisson 
equation for the pressure. The average pressure, (P^ + P^)/2, is 
computed at an interface vertex. From the average pressure and the pressure 
jump we can compute a pressure gradient centered on triangles, within and 


23 



without the droplet, for inclusion in the momentum equation. 

The radius of curvature is computed from a parametric cubic spline 
interpolant to the interface vertices. If we denote the interface 

vertices by r^=(x^, ) f i=1,...N, with r^=r^ , we define a pseudo arc length 

parameter, s, so that the spline knots occur at the points 

Si = ° 

s. = s. + Ir.-r. I i=2,...N. (18) 

We then generate the twice differentiable periodic spline interpolants x(s) 

and y(s) from the data s^^, x^, y^, i=1,...,N as prescribed by deBoor 

(1978). The curvature is then given by 

|x'y" - y’x"| 

(19) 

(x* ^ + y» 2) 3/2 

where a prime indicates differentiation with respect to the parameter s. The 
sign of R at an interface vertex, r^, is given by the sign of 

_ — — _ A A 

[(r. ,-r. ) X (r.-r. ,)].z, where z is the unit vector in the z direction. 

The parametric spline fit is also used for regridding. When the 
regridding algorithm calls for the bisection of a triangle side which borders 
the two media, a new vertex is added on the spline interpolant between the 
indicated vertices rather them bisecting the straight line segment. A 
straight line bisection introduces spurious interface oscillations (Foote, 
1973) whereas bisecting the spline maintains the general overall shape of the 
interface. 
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a. Droplet Oscillations 


In order to test our algorithm for surface tension, we performed 
calculations of droplets which oscillate under the effects of surface 
tension. A linear pheory for small amplitude oscillations on cylindrical 
jets was first given by Rayleigh (1879). When a perturbation is totally in 
the plane perpendicular to the axis of the cylinder, Rayleigh found the 
frequency, oi, for the oscillations is given by 

= (n3 - n) — ^ (20) 

pa^ 

where the surface of the jet is given in polar coordinates by 

r = a + ecos(n0). (21) 

Prom Equation (20) it is clear that the lowest oscillating mode is given by 

i 

n=2. Rayleigh used Equation (20) to interpret his experiments with jets. 

For large amplitude oscillations he found the experimental frequency to 
diverge from that predicated by his linear theory and attributed errors to 
nonlinear effects. 

In the numerical calculation presented we study an n = 2 oscillation. 

We have ta):en the parameters a = 0.0125 cm, and o = 30 dynes /cm, values which 
are typical for droplet combustion problems. We used a droplet density of 
2g/cm^ and an external fluid density of 1g/cm^. The results of a calculation 
with E = 0.2a = 0.0025 cm are shown in Figure 4. The numerical oscillation 
period is approximately 1.25x10”^ s. In order to compare this result with 
Rayleigh's theory, we must first correct his result for the effect of the 
presence of the external fluid. Equation (20) then becomes 
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a 


( 20 ' ) 


0^ = (n^-n)- r 3 

{ p +p )a-’ 
a e 


where p^ is the droplet density and p^ is the density of the external 
fluid. With the period de'fined cis 2ir/“» Equation (20') gives a period of 
1.13x10“^ s. The discrepancy between the numerical and theoretical results 
can be explained by the finite grid spacing. However, given Rayleigh's 
experience with large amplitude oscillations, it is reasonable to expect our 
computational period to differ somewhat from that given by the linear theory. 
Further calculations were performed with smaller amplitudes, e, to see if any 
of the difference is attributable to nonlinear perturbation effects or if the 
linear theory is directly applicable in this regime. In addition, different 
density ratios were used, viz. 10:1 and 800:1, which more closely 
approximated a combustion environment and which allowed the testing of the 
effects of the external fluid density on the numerical convergence of the 
pressure algorithm. The net result was that all the difference between 
theory and the numerical result is consistent with second-order convergence 
to the theoretical frequency for small perturbations and small grid size. 


b. Incompressible, Inviscid Flow about a Droplet with Surface Tension 

The second test of the surface tension algorithm was a recalculation of 
the initial benchmark problem, but with surface tension forces turned on. 
This test was necessary to check whether the code could allow for more 
radical interface deformations and whether the spline fit would properly 
allow the droplet to separate into smaller droplets or, alternatively, for 
many smaller droplets to coalesce. 
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Figure 5 shows the results of a calculation with surface tension for the 
same initial conditions as used in the calculation without surface tension 
(Figure 3), As in the case without surface tension, the internal droplet 
flow is driven by compression parallel to the external flow and is initially 
normal to the external flow. A recirculation zone is formed in the wake of 
the compressed droplet and the droplet deforms into a kidney shape between 
the opposing streams of the external flow and the recirculation zone. With 
the relatively large surface tension forces used in this calculation, further 
stretching of the droplet is curtailed. Instead of the droplet deforming 
into a film around the recirculation zone, the rear of the droplet begins to 
oscillate under the restoring force provided by surface tension, "nie 
oscillation arises at the rear of the droplet at a wavelength equal to the 
droplet diameter, TSie large deformations seen at later times have the 
shortest wavelengths which can be supported by the grid resolution at those 
times. These higher modes are excited numerically through wiggles induced by 
the spline fit to the interface vertices and by physical oscillations induced 
by the recirculating flow at the rear of the droplet. Ihe spline routines 
have been recoded for a higher order spline, but these algorithms have not 
been incorporated into the main routines at the time of this report. Ihe 
front of the droplet remains smooth throughout the calculation despite the 
large nonlinear oscillations occurring at the rear of the droplet, and the 
general droplet shape and behavior are consistent with experiments performed 
in low viscosity fluids. 

The Lagrangian pathlines for the vertices again show the development of 
a recirculation zone in the wake of the droplet. Initially this zone is 
similar to the one in the calculation without surface tension. The primary 
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effect of surface tension on the external flow is that oscillations are 
superimposed on the recirculating flow at the rear of the droplet. For the 
calculation shown here, this oscillation is sufficient to disrupt the regular 
flow pattern in the droplet wake and to induce higher mode oscillations. The 
effect of surface tension on the internal droplet flow appears in the 
retardation and cessation of droplet thinning around the recirculation zone 
and in the increased mixing due to the droplet oscillations. The internal 
flow remains laminar at the front of the droplet even in the presence of the 
large oscillations at the rear. The external and internal flow patterns and 
droplet shape at later times agree qualitatively with experimental shapes and 
flow patterns at high Reynolds number (Clift et. al. , 1978). This agreement 
extends to three dimensional droplets as well since experiments of bubbles 
and droplets between parallel plates show results similar to experiments of 
unconfined droplets and bubbles at their planes of symmetry. The calculation 
does not include viscosity so that the Reynolds number is large and limited 
only by an effective cell Reynolds number. 
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3. Viscosity 


The next step in the construction of the droplet combustion model is to 
include algorithms for viscosity and compressibility. The equations of 
motion of a viscous fluid in tvro dimensions require the additional terms 


dv 

— — = • • •+ V • vVv 
dt X 


dv 

— Z = • • •+ V • vVv 
dt y 


( 22 ) 


As discussed above, the formulation of the finite-difference algorithms 
required velocities to be specified at triangle centroids. Gradients of 

•> 

velocity components such as those in Equation (22) are therefore difficult to 
express to high order accuracy and the regions over which the approximations 
must be made are irregular and costly to compute. 

A numerical algorithm which is easier to implement can be derived by 
expressing the change in vorticity at a grid point due to viscosity as 

# = (23) 


since in two dimensions the vorticity is in the z direction. The easiest 
algorithm to implement is one which introduces the necessary changes in the 
triangle velocities about a vertex to satisfy Equation (23). If the choice 
is made to enforce equal contributions from each of the triangles about the 
vertex, then Equation (22) is satisfied. For example, consider an 
unperturbed shear flow parallel to the x-axis. For this flow Vv^^ defines 
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the shear profile and Vvy is zero. Therefore by Equation (22) only the x- 
components of velocity in regions adjacent to the shear layer should change 
and the y-components of velocity should remain unchanged. The choice of 
equal triangle contributions to the change in vorticity dictated by Equation 
(23) ensures that the numerical algorithm will induce change only in the x- 
velocities while the changes in y-velocities will be identically zero. In 
this sense, Equation (22) is used as a conservation law to ensure proper 
behavior of the finite difference algorithm. 

Although this algorithm was simpler to code, the specification of equal 
contributions from all triangles about a vertex was difficult to enforce 
except for regular grids. The determination of how the required changes in 
vorticity were to be translated to velocity changes was ambiguous for 
different grid geometries. Ihe algorithm produced the correct spreading 
rates for a shear profile, but only for very regular grid geometries. For an 
arbitrary grid a more detailed prescription was necessary. 

A discretization for which V^v is a triangle-centered quantity as in 
Equation (22) remains desirable. If in the finite-difference formulation for 
Equation (22) the coefficient of viscosity is centered on triangles, any 
ambiguity at interfaces is avoided for stratified fluids, whereas special 
algorithms would be needed for a vertex-centered coefficient of viscosity. 
TSiis placement of variables puts the viscosity on the same footing as the 
density. Temporal changes in the triangle velocities are straightforward to 
compute, since now 

!!t . 

dt \ 


( V2y) 


(24) 
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where the subscript "t" indicates that all quantities are triangle centered, 
a. Spreading of a Shear Layer. 

This algorithm was tested in a calculation of the spreading of a shear 
layer of initially zero thickness given by 

V (X, y, t=0) = + for y ^ y^ , 
where y^ is the original location of the vortex sheet. Ihe velocity 
distribution across this layer will evolve as 

V (x,y, t) = erf x (25) 

and the width Ay of the layer will grow as 

AY - 8( vt) ^ ''2^ 

For the test calculation the grid was initialized to center a vortex 
sheet in a grid 16 cells wide with an initial layer width of zero, llie two 
opposing streams had initially constant velocity profiles and the evolution 
of the interface between the streams was governed by the same algorithms as 
the interior of either fluid with no special interface boundary condition. 

The boundary condition at the sides of the computational region were 
periodic, and the top and bottom had free-slip boundary conditions. At the 
end of the calculation the layer width agreed exactly with the theory and the 
layer extended over the whole mesh. Ihe velocity profile for each stream 
coincided with that given by Equation (25) to within round-off error. Ihe y- 
components of the velocity remained zero indicating that the algorithm was 
working well for the grid distortions presented by the problem. 
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b. Incompressible Flow about a Droplet with Viscosity and Surface Tension. 

Before the droplet runs began two modifications were made in the code. 
The first was the addition of a new initial condition. All previous runs had 
used an impulsively started air flow, with the addition of viscosity, this 
led to large momentum transfers across the droplet interface early in the 
calculation. The new initial condition specifies a steady-state flow field 
derived from a streamf unction calculation. This provides a much smoother 
start and a closer representation of the actual physical conditions the 
droplet would see. The second modification was to the residual error 
algorithm which corrects for the effects of keeping straight triangle sides. 

A mistake was found which became evident only for large momentum transfers 
across an interface. The error was corrected and the problem was eliminated. 

The first viscous calculations were of air flow past a kerosene droplet 
and included the effects of both surface tension and viscosity. The physical 
parameters were appropriate for a combustor environment: 


density of kerosene 
density of air 
surface tension (STP) 
viscosity of kerosene 
viscosity of air 
air velocity 
droplet radius 


0.82 g/cm^ 

.0013 g/cm^ 

30 dynes /cm 
1.8 centipoise 
.018 centipoise 
100 m/sec 
125 microns. 


Figure 6 follows the evolution of the internal and external flow fields for 
the calculation. At an air velocity of 100 m/sec .and a droplet radius of 125 
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microns, the corresponding Reynolds number is roughly 1600. Boundary 
conditions for the computation are periodic on the sides of the computational 
region and reflective at the top and bottom. The passage of fluid through the 
mesh can be tracked by the pathlines of the uppermost and lowermost vertices 
next to the top and bottom of the computational region. These vertices are 
slightly above and below their neighbors due to the algorithm used to 
calculate the initial grid. Their position can be tracked through all nine 
frames, showing that the fluid originally behind the droplet has progessed 
through the mesh and has interacted with the face of the (next) droplet. Note 
the initial frame is not at t=0.0 in order to accumulate particle pathlines 
which are indicative of the originally laminar flov;. 

The first clear indication of the development of the recirculation region 
is seen in the fourth insert where a pair of counter-rotating vortices are 
evident. The recirculation zone continues to develop throughout the 
calculation, although at times the vortex pair is not as evident due to the, - 
deletion and addition of vertices which interrupts the continuity of the 
pathlines. By the last insert it appears that another pair of vortices isr 
forming near the droplet, indicating that the original pair may be shed. 

There is now large distortion of the leading face of the droplet, and the 
droplet is about to enter the wake of the preceding droplet. Distortions in 
the face of the droplet are evident by at least the seventh frame, and appear 
to be due to increased curvature and condensing of the streamlines in the 
external flow caused by the approaching wake. The internal velocities are 
small compared to the external flow rates and therefore cannot be 
distinguished as pathlines. However, an indication of the (small) internal 
recirculation may be obtain by comparing internal vertex positions at various 
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times teps. Figure 7 shows the grid at the corresponding times in the 
calculation. Vertex addition has evidently occurred primarily where needed, 
in the developing wake of the droplet and all along the droplet interface. 

Figure 8 shows the pathlines for a simulation with identical parameters 
except for the flow speed, which is increased to 120 m/sec for a Reynolds 
number of 2000. The fluid now completely passes through the mesh, with the 
fluid initially near the droplet completely passing the next droplet. 
Therefore the droplet has interacted with the weOce of the preceding droplet 
for one droplet diameter. The initial flow about the droplet is seen to be 
quite similar except for a more pronounced flattening of the face of the 
droplet due to the higher flow speed, 'ttie wake develops in much the same 
manner, but it now interacts strongly with the flow at the forward stagnation 
point on the droplet. Oscillations in the flow due to the wake are 
transmitted to the forward face of the droplet and give rise to fairly large 
perturbations. As seen in Figure 9, the computational grid is in need of 
further refinement at this time because the perturbations cannot be resolved 
by the length scales originally chosen for the run. One of the crests of the 
surface wave is gridded by a single triangle, a situation which allows no 
communication of that surface fluid with the interior of the droplet. In 
order to continue the simulation better resolution must be obtained about the 
droplet surface. A new algorithm has been developed to permit higher 
resolution near points of large curvature of material interfaces, but the 
algorithm was not implemented at the time of these calculations. 
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FUTURE DEVELOPMENTS 


The next step in the construction of the droplet combustion model is to 
include an algorithm for compressibility. The addition of compressibility 
will occur in two ways depending on the characteristic flew velocities -in the 
calculations. When flow speeds are slow with respect to the sound speed, we 
do not want the timestep to be limited by the Courant condition. In such 
cases the sound waves can he filtered out of the solution by altering the 
pressure field to account for local divergences on the time scale of the 
fluid flow (Jones and Boris, 1979). These divergences, which arise for 

i 

example, from heat release, are introduced into the pressure Poisson Equation 
in a manner similar to that for incompressible flov;. However, there is a re- 
striction that the relaxation times occur at the proper time scale. For the 
triangular mesh, such additions should be easy to implement since a 
divergence correction term is already used to account for the effects of 
maintaining straight triangle sides. 

In the case for which sound waves must be included, an energy evolution " 
equation and an equation of state must be included in the finite difference 
algorithms as well. The algorithm which will be used for the equation of 
state expresses the density as a function of the pressure and energy. Given 
a new internal energy derived from the energy evolution equation and an 
approximation to the pressure, density is calculated from the equation of 
state. This equation of state density is compared to the density derived 
from the fluid dynamics and the difference is iterated to zero. This 
solution is the inverse of the usual algorithms for the equation of state 
which express pressure as a function of density and energy. The method has 
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been tested extensively for a one-dimensional restructuring mesh in the code 
ADINC (Boris, 1979; Fritts et al., 1981). The technique works well in one 
dimension and exhibits diminished finite difference error propagation due to 
the fact that numerical errors in pressure result in small density 
fluctuations. In the usual technique, small density errors can give rise to 
large pressure fluctuations, and hence a larger error propagation. 

An energy evolution equation, 

=_eV «v-V • (pv) + V • XVt, (26) 

dt 

will also be needed to account for the effects of thermal conductivity, 
represented by the last term in Eq. (26). If both energy and temperature are 
carried as vertex centered quantities, then there are no new techniques 
required for the first and third terms, since the finite-difference 
approximations for similar terms are well tested. The center term must be 
carried as an average, since pressures are vertex centered while velocities 
are triangle centered. The incorporation of reactions is rather 
straightforward and will follow previously tested techniques given by Oran 
and Boris ( 1 981 ) . 

The three-dimensional analogue of a triangular grid is a tetrahedral grid 
in which surfaces are tessalated by triangles. Although the addition of one 
more dimension introduces new complications in the reconnection algorithms, 
much of what was learned from the two-dimensional case carries over intact 
into three dimensions. Vertices can still be deleted by successive 
reconnections to isolate a vertex within a single tetrahedron. At that 
point, the vertex and its four lines can be eliminated from the grid. 

Vertices cein be added within tetrahedra, in the plane of a triangle and on 
lines without major modification of the techniques used in two dimensions. 
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The conservation criteria used for reconnecting, adding and deleting cells in 
three dimensions usually involves either extending integrals to one higher 
dimension, or measuring an angle between planes rather than lines. 

Similarly, the hydrodynamics finite-difference algorithms are logical 
extensions of the two-dimensional algorithms. The use of primitive variables 
allows a simple extension for the vorticity integrals, and the solution of 
Poisson's Dguation still requires just one matrix inversion. 
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CONCLUSION 


In this report we have described the basic algorithms in the two- 
dimensional Lagrangian, incompressible Cartesian code, SPLISH. A main 
advantage of the Lagrangian technique is the general property of the 
restructuring triangular mesh, which allows reconnection of vertices without 
adding numerical diffusion. This technique is accurate at material 
interfaces even though the interfaces undergo convolutions and may evolve 
into multi -connected surfaces. Because of the potential advantages of such a 
technique to combustion problems, we have begun the process of converting the 
code for the study of flows in and around burning fuel droplets. 

We have described and presented benchmarks for two new algorithms which 
have been added to the basic fluid dynamics code: one for surface tension and 
one for viscosity. Surface tension is included as a ^ump in pressure across 
an interface by casting the surface tension forces in the form of a gradient 
of a potential. The algorithm has been benchmarked by comparing numerical 
solutions of the oscillations of an n = 2 normal mode to the results of an 
analytic solution. The difference between the exact and numerical solution 
becomes smaller as grid resolution is improved. The viscosity algorithm is 
presented and tested by calculating the spreading of a viscous shear layer 
and comparing this to the analytic solution. Here the analytic and numerical 
solutions are virtually identical. 

Finally, we have discussed initial calculations of the behavior of 
dense fuel droplets in a flowing gas. Droplet flows with and without surface 
tension and with and without viscosity are discussed. Calculations of 
kerosene droplets in air are presented. These show both internal and 
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external droplet flows as well as distortion of the droplet due to the 
relative flow. Also, we see how vortex pairs develop and are shed behind the 
droplet. Droplet-droplet interactions occur when the distorted flow induced 
from one droplet reaches another. Ihe results of the calculations are 
illustrated by sequences of frames from a computer generated movie of fluid 
particle positions. 

The restructuring mesh has been shown to be capable of accurately 
tracking interfaces despite transition to multiply connected flows. As a 
result calculations of distorting and shattering droplets can now be 
performed entirely from first principles, without recourse to approximations 
or phenomenological models. The numerical technique is appropriate to the 
tracking of flame fronts as well. There are a number of future directions 
that can be taken in the development of the model. Algorithms which make the 
code compressible have been developed and must be implemented. When this is 
done, we can consider effects such as thermal conduction and chemical 
reactions. Implementation of these algorithms in the cylindrical version of 
this code instead of the currently used cartesian version would allow the 
study of a round instead of cylindrical droplets. Some of the basic 
algorithms need to extend the code to three-dimensions have been worked out 
with tetrahedra replacing triangles. However, these need considerable 
development before they can be used here. 
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Fig. 2. Portions of a grid illustrating possible reconnections. In part (a) 
the dashed diagonal will be chosen for the shaded quadrilateral rather than 
the present, longer, diagonal. In part (b) the diagonal cannot be 
reconnected since the alternative diagonal, though shorter, lies outside the 
"inverted" quadrilateral. 
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Fig. 3 Frames from a computer generated movie of a simulation without 
surface tension or viscosity showino droplet deformation and shattering due 
to the external flow. 


TIME = 0.0 

3.00 X 10-< 

6.50x10-< 

MmiaTOMMMteaggjyiM 

snaraaaaXmmmwAVArjaa 

iTMBraaafi^^ 

kTaT^ATATA^ICtI^ATaT^^^^ 

unR0urATj9w^^^ 

■UfATATATAT^W^T^r^^AW^ 

■k^^TATATA^Cur^ATATSfU^^ 

■'iiTATSkTATATAC^A^^^^ 

■ k^TAigikTATAT4»w^^^ 
■'ATATATATAT^>'^I€^4^Ti0r^^^ 

■'iiTATinirA^^il^K'^^ 

■ kTATATATAW^Jv^^ 

■A^Ti»kTAnV^^ 

■ WATATATAlA^fVf!S8^^^ 
■kWATATiK^ICTl^ATA'I^^ 
BwiaaMkTATia^Tlw^ 

naasanBuir^-ff^M 

vjnirjBr^^^ 

amkiak^^^AlDl^iW^ 

kiakTA^^WisS 

»akiakn^Tv^<^>^WAT^^ 

kTAWATATA^^;^^^^ 

iTififAlSlkTATA^TAn^^ 

CrATATATA^^.^:^^>A^^^ 

*0S0S0ai!urA^^ 

■ UkTATATATAT^I;(ri^T^ 

■ kTAWATA^I^^^T^^TATA^TATi 
■'A^TATATiSk^Anv^ 

■ :raBwrasra)iii!>iifieieiaiiBra!asri 

■ kTATAlkTAW^;S»^^^^ 

■ ATATAnTATATi»ATATAt»V^ 
■rankTATak^^)^^^ 
■>ATiMakTiaA^^;y^^ 

BWi!IATAWJ8»!W^^ 

nirjiBxmir^AV^^ 

kTATATATA^A^CvS^^^^^ 

vaxaaa^^ 

nsrjiairM^^^ 

'aTatataWa^^^^ 

yA^AWg^gjgl^”^TOkTATO 

'Meuiegfak^Ar^^^ 

klUIATAiakWj^W 

ATATiR^TAT^X^I^I^^AWAWA 

kTATATATATA^WA^VA^^^^^ 

niTATATATAri^^I^A^ia^^ 

iTATATS0!fi^^^liB^^^ 

MAWBi«gglBli?5»«ii!riW 

wmy^CTiS&llW&gwi^witw 

MMSOATAlS^TfT*^ 

USkTAWATiBtakTAT^^^ 

■ aiTATA^TAT^Xtlt^^^ 

■ kTiOkTATA^I^^X^^ATAiakT^^^ 
MnniimsrA^^^^ 
■kTATakT0.fir^%M^€tA^ 

m^MsaBis9c<\mmti^iawsixa 

■anK0ATjWl»AWA^^ 

IkTA^TiirAW^W^SBi^^ 

■wakiATaktaA^fffi^^ 

lak^^tATATAt^^TATj^T^^ 

wjomin^AYA'rM 

VrjSUFOOk&A^AT^ 

neum^ArAW^^^ 

ATAis»i;ve'AB^^T^ai9ii^^^^ 

ATATifla^\?^|!%V^^ 

0uriiiAWAVi»ATAnTi^^^^ 

’rasaBOBOBll^ 

0!ff£nsrM»MAyjBavjsr£^^ 


Fig. 4 Computer generated frames from a Lagrangian simulation of an 
oscillating droplet for two periods of an g = 2 normal mode. 
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Fig. 5 Frames from a simulation with surface tension but no viscosity 
showing droplet deformation due to the external flow. 


Fig. 6 Pathlines of the fluid flow from a computer generated movie of 
incompressible air flow past a deforming kerosene droplet. Heads of the 
pathlines are the current vertex positions and the tails are made up of the 
previous six positions. Ihe flow speed is 100 m/sec. and the effects of 
surface tension and viscosity are included. 
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Fig. 7 Frames showing the triangular grid at the same times as shown for 
the pathlines in Figure 6. 
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Fig. 8 Pathlines for the fluid flew when the flow speed is 120 m/sec. All 
other parameters for the calculation are identical with those for the 
simulation shown in Figure 6. 
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Fig. 9 The restructuring triangular grid at the times corresponding to the 
frames in Figure 8. 
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